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The properties of the reduced density matrix describing an interval of N sites in 
an infinite chain of free electrons are investigated. A commuting operator is found for 
arbitrary filling and also for open chains. For a half filled periodic chain it is used to 
determine the eigenfunctions for the dominant eigenvalues analytically in the continuum 
limit. Relations to the critical six-vertex model are discussed. 



Reduced density matrices for some portion of a larger system play a central role in the density-matrix 
renormalization group (DMRG) method [1-3]. For this reason, they have been investigated for a number 
of model systems in the last years [4-8] . It was found that for free electrons or bosons the reduced density 
matrices have exponential form exp(—H.) and thus look like thermal operators. It was also realized that 
one can base the determination of Ti on the one-particle correlation functions in the state one is study- 
ing [7,9-11]. However, for critical lattice systems, which are particularly interesting in view of DMRG 
applications, the properties of TL have only been studied numerically so far. For free electrons hopping 
on a chain one can see in this way the typical finite-size effects on the spectra and the concentration of 
the eigenfunctions near the boundaries. Nevertheless, a more explicit solution, which also makes contact 
with some results found in field theory [12] in calculations of entanglement entropies, would be desirable. 

In the following we present such a solution by looking at the hopping model in a somewhat different 
way. We show that there exists a relatively simple operator T which commutes with the reduced density 
matrix for the ground state and thus has the same eigenfunctions. This holds for arbitrary filling and 
also for a subsystem at the end of a semi-infinite chain. The case of a half-filled periodic chain is then 
studied in detail. Working with T and taking a proper continuum limit allows to determine the general 
character of the single-particle eigenstates. For the low-lying states it is found that the eigenvalues of 
Ti and T even coincide up to a scale factor. We also discuss the connection with the critical six-vertex 
model. This allows to derive the spectrum of Ti by a conformal mapping. It also allows to view the 
density matrix as a particular transfer matrix and the commutation relation as a special case of similar 
relations which occur in the treatment of integrable two-dimensional models. 



We consider a system of free fermions hopping between neighbouring sites of an infinite linear chain. 
The corresponding Hamiltonian reads 



where t is the hopping matrix element and the 'hat' denotes quantities of the total system. The ground 
state |0 > is a Slater determinant corresponding to a certain filling factor n. As subsystem we take N 
consecutive sites, labelled by i,j = 1, 2, ..N. The reduced density matrix p has the general form 



I. BASIC FORMULAE 
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p = JCcxp (-H) 



where K, is a normalization constant and 
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The matrix is completely determined by the one-particle correlation function of the total system 

C mn =< 0| 4c |0 > (4) 
which in our case is given by C mn — C(m — n) where 

(5) 



C{m) = f ^ n(q) e^ m 

and n(q) = 1 for \q\ < ttu and zero otherwise. This leads to 

„, . sinfirnm) 
C(m) = 
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If C is the N x N submatrix of C where the sites are restricted to the subsystem, one has [7,9] 

H = ln [(1-C)/C] (7) 

Thus H and C have common eigenfunctions and their eigenvalues Ek and (k are related according to 
(7). The (k lie between and 1 while the vary between — oo and oo. For large N, the spectrum of 
C must approach that of C which is given by the step function n{q). For this reason, most of the (k lie 
exponentially close to and 1 [11,8] and C is not easy to handle numerically. 



II. COMMUTING MATRICES 

This difficulty can be circumvented by working with a matrix which commutes with C (and thus H), 
but has a simpler spectrum and also a simpler form. The existence of such a matrix is suggested by the 
structure of Hij obtained from numerical calculations. One finds that Hij has approximately, but not 
exactly, tridiagonal form. A similar situation occurs for the transfer matrix of the two-dimensional Ising 
model, and there a commuting tridiagonal matrix exists [13]. Trying therefore 



( d\ ti 

t\ d 2 t<i 
h d 3 



\ 



*3 



(8) 



one finds indeed that [C, T] = 0, if the coefficients are chosen as 
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The mirror symmetry of the subsystem is reflected directly in these coefficients. Their linear increase 
near the boundaries is analogous to the situation found for corner transfer matrices [14-17]. Physically, 
the matrix T describes a hopping model with Hamiltonian 
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where both the hopping and the single-site energies increase towards the middle of the subsystem. In the 
case of half filling (n = 1/2), the diagonal terms are zero and one has only hopping. It should be noted 
that T is only determined up to a constant factor, and we have chosen a normalization where the largest 
coefficients are of order one. 

A similar result holds if one considers a total system which is semi-infinite and chooses as subsystem 
the first N sites next to the end. This geometry corresponds to the one in DMRG calculations with open 
boundaries. In this case, one has to replace Cij = C(i — j) by the expression C-j = C(i — j) — C(i + j) 
in all formulae. Then one finds for half filling, that [C S ,T S ] — if the ti are given by 
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(11) 



The hopping varies again linearly near the interface with the rest of the system and is maximal at the 
open boundary. Basically, T s corresponds to the right half of the system described by T which is easy to 
understand in terms of the geometry. In the following we study T for the case of a half-filled system. 



III. EIGENVECTORS 



We now use T to determine the common eigenvectors of C, H and T . From numerical calculations 
one sees that for the eigenvalues A of T which are smallest in magnitude the eigenvectors ip show rapid 
oscillations. It is convenient to take these out and also to separate the amplitudes at even and odd sites 
by writing 

= (-iy +1 V2j -u = (-i) J+ V 2j (12) 

This leads to the coupled equations 

-t 2j -2^j-i + *2i-i*i = A*,, tij-i&j - tijQj+i - A*,, (13) 

If N is odd, one eigenstate can be found exactly. Then A = is an eigenvalue and the equations decouple. 
One finds that ^fj — 0, while $j follows from the the recursion relation 

$ J+1 = ^I$„ (14) 
t 2 j 

which can be solved explicitly. This " central state" (the A are distributed symmetrically with respect to 
zero) follows from the structure of the matrix and exists also for T s . It has also been found in the study 
of Harper's equation where a similar matrix (with sinusoidal elements) appears [19]. To see the form 
of one can consider large N and take a continuum limit by writing (2j — 1)/N — 1/2N = x, t 2 j-i = 
t(x),$j — $(x). Expanding the quantities, one then finds the differential equation 



with the solution 



VHx) \/x(l — x) 
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which reflects directly the form of the hopping amplitude t(x). This state was studied numerically in [8]. 
The same continuum limit can also be used to determine the low-lying eigenstates in general. The 
equations (13) then become 

(2t±+t>Wx)=^(x), -(2t±+t'Mx)=^(x), (17) 
where a = XN absorbs the 1/N dependence of the tj. Combining the equations and putting 

$(a;) = -^=L (18) 



one finds for \ 

dx dx 
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^-r-t—x = n\ (19) 



The substitution 

a — . , , 

2 y l-x h 



«=^»(t £ t). ( 2 °) 



then leads to the oscillator equation for x{u 

i. 

du 2 

The general form of $ is therefore 



S+M 2 X = 0. (21) 



$(x) = — sin 

^x(l -x) 



l ln{ r~x )+a 



(22) 



The function ^f(x) is obtained by replacing the sine with minus cosine according to the relations (17). 
Such logarithmic oscillations of the eigenfunctions were also noted in [20] and had been found earlier for 
corner transfer matrices [21-23] where matrices with linearly increasing elements are involved. For the 
Ising model, the square root appears there, too, and reflects the surface exponent x s = 1/2, while for 
the Gaussian model there is no such prefactor. To determine the allowed values for /i, one needs proper 
boundary conditions. From (13) one has tx^i = A$i and ti^N/ 2 = These lead to the values 

a = ±7r/4 for the phase and asymptotically one finds 

^ = ±^£^(2^ + 1), k = 0,1,2,.... (23) 

i.e. equidistant levels scaling as 1/lnN. For a bosonic field, such a spectrum was also found in [20]. 
However, this behaviour can be seen only for very large N. For moderate N of the order 10-100 or even 
1000, the dispersion relation Uk vs. k still shows some curvature [6,8] and the eigenvalues are found 
to vary as l/(lnN + bk) rather than as 1/lnN. The constant bo for the lowest state, for example, is 
approximately 2.6. Therefore one has considerable finite-size corrections to the asymptotic law. These 
features are also known from studies of corner transfer matrices [18,21]. 

The trigonometric functions % vanish near one end of the subsystem and are maximal near the other 
one. For each eigenvalue, they are related by symmetry such that $(x) = ±*(1 — x). 
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IV. DIRECT TREATMENT OF C 



With the experience from the treatment of T one can also obtain the eigenvectors directly from C. 
Since C is not a sparse matrix, this leads to an integral equation. Working again with $ and 4", one has 
instead of (13) the equations 

J2 D(2* - 2.7 + = C *i, - E D ( 21 - 2 -1 - = C *i> (24) 

where D(l) = l/(ir I) is the denominator of C(l) and ( = ( — C(0) = £ — 1/2. Combining these, gives a 
single equation of the form 

E^^=C 2 $i (25) 

Before taking the continuum limit, one has to single out the diagonal terms Kn which are positive, while 
the others are negative. Since 

1 N/2 1 

^ = ^E [20 - i) + 1]a ( 26 ) 

depends only weakly on i and N, one approximates it by its limit for large i and N, Ku = K = 1/4. The 
equation in the continuum limit then becomes 

r 1 - i 

/ dx'K(x, x 1 ) ${x') = (C 2 - -) (27) 



where 



X(a;,a;') 



2tt 2 a; - 

Using the same substitutions as before, one obtains 



In— Zn- 



1 — .t 1 — x 1 



(28) 



2 71 " J-oo sti(u — u) 4 

Since one has a difference kernel now, the equation can be solved by Fourier transformation. Writing \ 
as a trigonometric function of argument jiu as in the solution of (21), one finds for ( 



C-5 



l + th{^) (30) 
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The eigenvalues e of H then follow from e = ln((l — C)/C) > which gives the simple result 

e = nu (31) 

i.e. up to a factor of n they are the same as the eigenvalues of NT. For the low-lying states and in the 
continuum limit, therefore, the matrices H and T are proportional to each other 

H = ttNT (32) 

This can be checked numerically, and one finds indeed, that e and A fulfill this relation with high accuracy. 
For TV = 16, for example, the relative deviations are of the order of 10~ 3 for the lowest states. In 
general, however, the relation (32) holds only approximately, and the deviations increase for the larger 
eigenvalues. The close relation of the two matrices can be understood from the structure of H which, to 
a first approximation, has the same bidiagonal form as T with very similar matrix elements. 
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V. CONNECTION TO TWO DIMENSIONS 



Using the Jordan- Wigner transformation, TC can be expressed in terms of spin one-half operators and 
becomes the Hamiltonian of an XX chain 

H = -2t + <Tlal +1 ) (33) 

n 

where a*, are Pauli matrices. It is known that this operator commutes with the row transfer matrix V 
of the six- vertex model with A = [24,14]. This model also corresponds to two uncoupled Ising lattices at 
their critical temperature. Moreover, the ground state of H is the state with maximal eigenvalue for the 
transfer matrix. It follows that one can interpret the reduced density matrix p as the partition function 
of such a system, having infinite size and a cut in it, along which the spins are kept as free variables 
[4,25]. In previous work on non-critical models, this cut was considered as half- infinite and one was lead 
directly to standard corner transfer matrices [4,5]. In the present case, the cut has a finite length N. 



This property of p can be used to redcrivc the low-lying spectrum of TL. For this purpose, one considers 
the Ricmann plane which is obtained by stacking such systems and joining different sides of the cut. If 
the cut is between x = — a and x = a, this is the Riemann plane of the function ln(z — a)/(z + a). 
For correlations between points on the cut, but in different sheets, p plays the role of a transfer matrix, 
and the asymptotic behaviour of such correlation functions is determined by the largest eigenvalues of p. 
Since one is at criticality, one can use the conformal transformation 

■ z — a . 

w = In (34) 

z + a y ' 

to map the complete Riemann plane to the full w plane. If one further cuts out small circles with radius 
b << a around the endpoints z — ±a of the cut, the map is to a vertical strip — uq < u < uo in the 
w plane, where w = u + iv and uq — ln(2a/b). A vertical correlation function in the strip between 
points lying Av = 2tt apart correponds to correlations between equivalent points in consecutive Riemann 
sheets and involves one power of p there. This allows to determine the spectrum of p (or H.) from the 
correlations in the strip as in [18]. For free boundary conditions along the small circles, this gives the 
positive e as 

£k = 2TTy(x s + k), k = 0,1,2... (35) 

where L — 2u is the width of the strip. Inserting the value x s = 1/2 for the surface exponent of the 
six-vertex model one then has 

e k = Jl £- ir){ 2k + l) t k = 0,1,2... (36) 

If one identifies 2a with TV and sets 6=1, this is the continuum result (23), (31) found previously. By 
choosing b < 1, one could also obtain an additive correction to InN. 

One should mention that a similar mapping was used in [26] to calculate entanglement entropies in 
the context of black holes. In this case, one needs p to calculate the entropy via S — —tr(p Inp). The 
result was that asymptotically S diverges as S — (c/3) InN, where c (equal to one here) is the conformal 
anomaly. This behaviour was verified recently for various spin chains [10,11,27] and can be understood 
in terms of the logarithmic dependence of Sk on the size N [20] . A discretized spectrum as in (36) was 
found previously for finite-size corner transfer matrices. Then 2a corresponds to the outer radius of the 
system and b to the inner one [18]. 
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The same approach can also be used to discuss p for the case of a subsystem at the end of a half-infinite 
chain. If the end is at x = 0, one is lead to a half-plane x > with a cut from x = to x — a, i.e. from the 
left edge into the system. This maps to the left half of the previous strip under (34). Therefore, one has 
to multiply et in (36) by a factor of 2. Such an increase of the single-particle eigenvalues (corresponding 
to a faster drop of the spectrum of p itself) can also be seen in numerical calculations on the lattice 
and is responsible for the better performance of the DMRG with open boundaries in the present case 
of a critical system. In the entanglement entropy, this gives a factor of 1/2 which one can interpret as 
resulting from the reduced number of interfaces between the subsystem and the rest. 



VI. DISCUSSION 



We have considered the reduced density matrix for a chain of free electrons in its ground state and have 
found results in two directions. On the one hand, we could derive the structure of the low-lying eigen- 
functions which are most important in DMRG calculations and also in general. These have their largest 
amplitudes near the interface(s) between the subsystem and its surrounding. Away from the interface, 
they decay with a power law involving the surface exponent x s = 1/2 of the six- vertex model. Thus 
the influence of the environment propagates far into the interior, as expected for a critical system. The 
spectrum of H was also obtained and checked by a conformal mapping in the associated two-dimensional 
model. These results complement previous numerical studies [6,8]. 

On a more general level we found the operator T, Eqn. (10), which commutes with the reduced density 
matrix p. We used it here as a tool for obtaining the eigenfunctions of p, but it is interesting in itself. 
Written in terms of Pauli matrices, it has a form analogous to (33), namely (for half filling) 

JV-l 

T = -2 5>(<^i+"Mn) (37) 

i=l 

and describes a spatially inhomogeneous XX chain of finite length. The commutation relation [p, T] = 
is completely analogous to the relation [V, H] = between the row transfer matrix V of the six-vertex 
model and the Hamiltonian H. of the total chain. It would be interesting, and give further insight, to 
derive it in a different way. One should mention that reduced density matrices have actually been used 
in studies of integrable models via vertex operators or via the algebraic Bethe ansatz. In particular, the 
six- vertex model and the equivalent XXZ spin chain were treated in [28,29]. This lead to exact, but rather 
complicated expressions for correlation functions along the cut. By contrast, the commuting operator 
gives a relatively clear physical picture of the problem under consideration. 
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